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1. INTRODUCTION 




The space shuttle main engine (SSME) has extremely complex internal flow structure. 
The geometry of the flow domain is three-dimensional with complicated topology. The flow 
is compressible, viscous and turbulent with large gradients in flow quantities and regions of 
recirculations. In recent years computer codes are being developed^ to solve the flow 
equations in different regions of the SSME such as the hot gas manifold (HGM) region. 
The analysis of the flow field in SSME involves several tedious steps. One is the geometrical 
modelling of the particular zone of the SSME being studied. It is usually available in the 
form of engineering drawings, in terms of algebraic equations for different pieces of the 
surfaces or in a CAD (computer aided design) system. Accessing the geometry definition, 
digitizing it and developing surface interpolations suitable for an interior grid generator 
requires considerable amount of manual effort. There are several types of grid generators 
available with some general-purpose finite element programs, such as NASTRAN, ADINA, 
ABQUS, etc. However, these programs require considerable amount of effort on the part 
of the user to input the geometry to the grid generators; also, the grid generated by those 
programs are not always the most appropriate grids for the flows being modelled. Next, an 
efficient and robust computational scheme for solving 3D Navier-Strokes equations has to 
be implemented for this class of problems. Post processing software has to be adapted to 
visualize and analyze the computed 3D flow field. Different elements of the above process 
have been studied in the past and other parts are yet to be developed. The current report 
discusses the progress made in a project to develop software for the analysis of the flow in 
the space shuttle main engine and similar complex internal flows. 

A CFD code for practical applications should have the following features. It should 
be reasonably accurate for the cleiss of problems it is designed to solve, with grids that 
can be accommodated on the present day computers. It should be robust in the sense 
that it is numerically stable for a broad range of initial and boundary conditions and 
geometrical parameters, and tolerate some variations in the grid resolution and structure. 
It should be computationally efficient for obtaining accurate solutions with reasonable 
computational and human resources. Standard of efficiency, however, is relative and it 
can only be measured against the current CFD software or which can be foreseen in the 
immediate future. Another important aspect of a CFD code is its usability, as to how 
much effort a user has to expend to solve practical problems with it. 
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For computing the viscous compressible flow inside the main engine where the flow 
undergoes complex turns through various chambers and ducts, it is necessary to discretize 
the physical space with several competing requirements. The geometry of the internal 
surfaces is typically represented in a CAD system or in some equivalent form by the 
designer. The surface data representation should be interfaced with suitable interpolation 
software. The refined spline surface representation of the flow boundaries will be the input 
for the grid generation routines. The topology of the grid structure depends on the flow 
solver algorithm to be used. Finite difference codes usually impose constraints on the grid 
structure such as the separability of the indices for efficient computational procedures, 
while the finite element method can be implemented with less stringent requirements on 
the grid structure. The grid should provide reasonable resolution of the flow field within 
the limits of the grid selected by the user. This requires providing more grid points and/or 
special methods in regions of large gradients of flow quantities, such as the viscous zones 
near solid boundaries. The gird should meet certain smoothness requirements so that the 
metrics of the curvilinear grid can be computed numerically and the computed metrics 
.are nonsingular. Unreasonably skewed grid cells or elements, and singular points in the 
grid where the local transformation of the physical space to computational space has very 
small or very large Jacobians, should be avoided if at all possible. Otherwise such grids will 
require special handling by the flow solver algorithm and also may give rise to numerical 
inaccuracies and instabilities. 

There are several grid generation techniques and special purpose codes which can gen- 
erate reasonable grids for simple two-dimensional and three-dimensional geometries, for 
both internal and external flows. These techniques fall under two classes: algebraic gener- 
ators and elliptic generators. Algebraic generators use various interpolation and stretching 
functions while elliptic generators solve a set of elliptic partial differential equations. While 
both techniques are effective for simple geometric regions, it is usually difficult to use them 
to develop a composite grid over a complex internal flow domain. Finite element commu- 
nity have developed extensive amount of software for generating algebraic grid suitable for 
finite element solvers. For example, NASTRAN (a general purpose finite element program 
primarily developed for structural analysis) contains grid generators for 2D and 3D struc- 
tures. Also, the program PATRAN (developed by PDA Engineering) contains 2D and 3D 
grid generators and pre- and post processing capabilities. In the current project some 
parts of the software such as PATRAN will be adapted and developed to generate body 
conforming, curvilinear finite element meshes of the flow domains inside the SSME. 
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Computation of the flow field inside the space shuttle main engine requires the ap- 
plication of the state-of-the-art CFD technology. Several computer codes^^"'*^ are under 
development to solve three dimensional Navier-Stokes equations with different turbulence 
models for analyzing the SSME internal flow, such as the flow through the how gas man- 
ifold (HGM). The computational methods^®”®^ used in the Navier-Stokes codes fall into 
two major categories: finite difference and finite element methods. Some of the algo- 
rithms are designed to solve th« unsteady compressible Navier-Stokes equations, either 
by explicit or by implicit factorization methods, using several hundred or thousands of 
time steps to reach, a steady-state solution cisymptotically. Other algorithms attempt to 
solve the steady-state equations by relaxation methods. All of them require body-fitting 
curvilinear grids with sufficient resolution. Grid requirements, however, differ greatly with 
the region being modelled and the algorithm used. Implicit factorization based on finite 
differences typically use global numerical transformations whereby the transformed grid in 
the computational space is uniform and rectilinear. This requires the grid to have indices 
which axe separable in the three directions for three dimensional problems, and also be 
reasonably smooth. However, such requirements may introduce grid singularities when 
complicated domains are discretized. Flow solver algorithm will have to deal with such 
grid singularities. Explicit schemes and finite element algorithms have less stringent re- 
quirements on the grid structure. However, explicit schemes are slow to converge because 
of the stability limitations on time step, particularly for large scale viscous problems. 

The finite element method is characterized by three basic features which are credited 
for the enormous success, the method ha.s enjoyed in the solution of practical engineering 
problems^®^ The first feature is that every computational domain is viewed as a collection 
of simple subdomains, called finite elements. This feature allows us to represent compli- 
cated geometries as assemblages of simple parts. It is a desirable feature in the solution of 
flow problems in complex configurations, not only to describe the complex geometry but 
also to choose the most suitable computational grid for a particular flow. This feature also 
allows us to place or remove any obstructions routinely into the flow field. The second fea- 
ture is that over each element the solution is represented by polynomials of desired degree. 
This allows us to compute the solution as a continuous function of position instead of at 
selected few points. Desired degree of approximation (e.g., linear, quadratic, etc.) can be 
easily and routinely specified without rewriting the whole or parts of the program. The 
third feature is that the relationship (i.e., the algebraic equations) between the solution 
and its dual variables (i.e., velocities and forces) is developed using a variational method, 
such as the Galerkin method. The boundary conditions are then applied on the algebraic 
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equations directly before solving. The three features of the finite element method also 
allow the easy development and interfacing of pre- and post-processors, and user-defined 
subroutines for equations for state and turbulence models. 

The Galerkin finite element method (i.e., the weight functions are the same as the 
approximation functions) applied to fiow problems always results in implicit schemes. The 
weighted-residual (or Petrov-Galerkin) method, in which the weight functions are differ- 
ent from the approximation functions, can be used in conjunction with explicit schemes to 
obtain explicit final equations. For example, by selecting the weight functions to be orthog- 
onal to the approximation functions, the mass matrix can be diagonalized. However, such 
considerations are entirely in the interest of obtaining explicit schemes and not necessarily 
in the interest of accuracy or even computational efficiency. In the current project implicit 
finite element scheme with suitable dissipation terms for stability is being developed. A 
relaxation procedure, known as the locally implicit scheme is being developed to solve the 
coupled set of algebraic equations efficiently. 

In the following sections we discuss the technical approach to the development of the 
finite element scheme and the relaxation procedure. Appendix I contains the details of the 
equations derived and Appendix II has a listing of the three dimensional finite element code 
for the compressible Navier-Stokes equations. Future reports will discuss the numerical 
results for specific problems. 


2. TECHNICAL APPROACH 
2.1 GOVERNING EQUATIONS 

In an Eulerian description, used most extensively in fluid dynamics, the coordinate 
system is fixed in space rather than in the body, and measurements of density, velocity, 
pressure, etc. are made for the material particle that happens to be in a given location 
at that particular time. The basic equations of a continuous medium in the Eulerian 
description are; 

Continuity Equation . - The law of conservation of mass leads to 

( 1 ) 

where p is the density of the medium, v is the velocity vector and x = (xi,X 2 ,X 3 ) the 
spatial coordinates. 
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Equations of Motion . - The law of balance of linear momentum leads to the celebrated 
Eulerian equation of motion, 


— (pv) + V • ipvv) = V 


. i Z? 


( 2 ) 


Here F the body force vector (measured per unit volume) and a the total stress tensor, 
which caji be divided into hydrostatic and viscous parts: 

2. = -pL + L (3) 


Here p denotes the hydrostatic pressure, r the viscous (or shear) stress tensor, and I 
denotes the unit tensor. 

An application of the law of balance of angular momentum and neglect of microstruc- 
tural effects such as couple stresses lead to the symmetry of stress tensor, 

<^ij = <^ji, Tij = Tji (a = ^) (4) 


Enersv Equation . - The law of conservation of energy (the first law of thermodynam- 
ics) leads to 

— (pe) + V • (pev) = V • (a ■ v) -t- F ■ u -t- p5 - V • ^ (5) 

at “ 

where e is the total energy per unit ma^s, 

1 

e = a ^ V • V 

2 ~ ~ 

£ being the specific internal energy, S is the rate of internal heat generation per unit mass, 
eind q is the heat flux vector or the rate of heat flow per unit area across the surface in the 
direction of its unit outward normal. 

Constitutive Equations . - The thermodynamic pressure p is related to the specific 
internal energy s and the density p through an equation of state, 

p = p{s,p) ( 6 ) 

and the viscous stress is related to the deformation rate tensor d through a constitutive 
equation of the form 



where 


+ ( 8 ) 

and c is the tensor of viscosities. 

For isotropic fluids obeying linear stress-strain relations (i.e. Newtonian fluids) we 

have 

T{j = 2/idt'y (9) 

where fx, is the viscosity. 

Initial Conditions . - At time t = 0, values of all the dependent variables {p,v,e,p) 
must be specified in the entire domain. It is not essential to specify all of these quantities 
at the same set of points. 

Boundary Conditions . - Depending on the type of the boundary (e.g., rigid boundary, 
free surface, interface, plane of symmetry, etc.), there are different kinds of boundary 

conditions in a problem. At a rigid .boundary, the normal component of the particle velocity 

must coincide with the normal component of the velocity of the rigid boundary. For a fixed 
(in time) boundary, the normal component of the particle velocity must be zero at that 
boundary. A plane of symmetry can be interpreted as a fixed boundary. On a free surface, 
the pressure must vanish. At an interface (and at a contact discontinuity) the pressure and 
the normal component of particle velocity must be continuous, and the density, internal 
energy and the tangential component of particle velocity may be discontinuous (i.e. jumps 
may occur). Across moving shock fronts, the Rankine-Hugoniot relations must be satisfied. 

2.2 FINITE ELExMENT MODEL 


Writing the governing equations in terms of the velocities, pressure, density and in- 
ternal energy, we obtain 

|-v.(«) = o 

^ (pu) -r V • (pvv) Vp = pV ■ d-r F 
at 

Q 

— (pe) -I- V • [pey) -r- V • (pu) = pV • {d- y) F- y — ^ ■ q (10) 

and p is given by the equation of state. If we assume that the body force, heat flux, and 
the internal heat generation are zero, the last two terms in the energy equation dropout. 
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For simplicity and computational convenience, we denote 


PU. = Y.i pe = E 

so that (10) become 

|W + V.K = 0 

^(£) + V • (vVl) + vP = AtV • d 
at 

-h V . {Ev) + V • (Pu) = mV . (d . u) 
at 

We seek approximate solutions to Eq. (11) using the finite element method. 


(11) 


Spatial Approximation . - Finite element approximations to Eq. (11) are sought over 
a typical element fl*: 

JV 

J = l 

N 

e = '£ (12) 

>=1 


where Vj{x) are the interpolation functions in space, pj, , and Ej are the unknown, 
time-dependent, nodal values to be determined. In Eq. (12) we have assumed for simplicity 
the same type (linear or quadratic) of interpolation functions for all the variables. The 
Galerkin approximation amounts to seeking solutions to Eq. (11) in the form (12) by 
making the errors in Eq. (11) orthogonal to the trial functions. This leads to the following 
local set of nonlinear ordinary differential equations in time. 


iA|W + [Bl{V>=0 
U|{£} + !M|{E> = {R} 


Here the superposed dot denote total differentiation with respect to time, and 

A{j — f ~ f a ’ 

Jq> Jq' dXk 


( 13 ) 
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Nij = [ {Vk^j) dx+ f 

•/n* jt=i •'^'’ 

r ^ ^ Q 

Mij = yz ^ — {Vki^j) dx 

Jo‘ fr; dxk 

Qik — f dx^ R-i - f ip{ ^ (p^fc) "I" [ /^'^4^id^dx (14) 

Jn« 7n* in* 


L 


where / denotes integration over the element volume. 


Equations (13) are to be further approximated (or numerically integrated with respect 
to time) to obtain a set of simultaneous algebraic equations. 


Temporal Approximations . - Equations (13) are of the general form 


(15) 


We approximate U(t) by 

n 

m=l,2,...,M (16) 

j=i 

where d>j (i) are approximation functions in time. Here we assume that <t>j are linear in t 
(i.e., n = 2): 

Mi) = ’ *^2(i) = ^, 0<t<At 

where At denotes the time increment. Then the time derivative of U is given by 

U={U2-Ui)/At (17) 


It cam be readily interpreted that Ui is the value of U at time t = n{At), and U 2 is the 
value o( U at t = {n l)At. Substituting Eq. (16) and (17) into Eq. (15), multiplying 
with ^ 2(0 and integrating over 0 to At. we obtain 

{t-'a} (18) 

Thus the unknown vector {Un->~i} can be solved in terms of the known vector ft 

should be noted that the temporal approximations (18) can be applied to the local set (13). 
There are other methods of time integration which can be incorporated into the code. 


~AtB 
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Equations (18) can be assembled in the usual manner to obtain the global equations, 
which must be solved iteratively (after imposing the initial and boundary conditions of the 
problem) for the nodal values, as the resulting algebraic equations are nonlinear. A flow 
chart of the computer program based on the formulation presented above is shown in Fig. 1. 
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Fig. 1, Flow Chart of the Computer Program 
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2.3 LOCALLY IMPLICIT APPROXIMATIONS 


For large problems, it is not possible to solve the global (linearized) equations by 
direct methods. An efficient iterative method of solution has been formulated and it is 
known as the locally implicit method^^^. This is based on a modified Gauss-Seidel iteration 
technique with, a symmetric inner iteration. 

The linearized equations (18) for an element j can be written in the form 

Lj AUj = Res (C/") + ^ Lfc + Q (19) 

where AC/j = — UJ and the summation on the right hand side of (19) is limited to 

the elements surrounding the element y, with which the finite element equations over the 
element j are coupled with. Equations (19) are solved by an iteration 

LMj 6 AUj = Res [U^] + ^ - Ly + Q (20) 

k^j 

where A(7”’+^ = ^6 AUj, LMj is a modification to the matrix Lj so as to achieve 

stability and rapid convergence of the iteration process. AU^} denotes either 

f )7l) 

or AUj^ ' depending on the latest available iterates for AC/*. The iteration process of the 
equation (20) is carried out starting at a different comer of the computational space for each 
iteration. Eight such iterations complete one symmetric modified Gauss-Seidel iteration 
per time step for 3-dimensional problems. This is a stable process with fast convergence 
properties in a local sense. It amounts to solving the equations (15) implicitly in a local 
sense for each node. It is not necessary to achieve full convergence at each time step if 
we need only the steady state solution. One symmetric sweep per time step is adequate. 
This process has been tested over a variety of model equations such as the 3-dimensional 
Poisson equation and one dimensional Burger’s equation. The same procedure has also 
been shown to work for two dimensional Euler equations with finite volume discretizations 
and artificial dissipation terms. 
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Appendix I 


Finite Element equations for Navier-Stokes Equations 

Variational formulation over an element for the Navier-Stokes equations in non- 
conservation form: 


0 = 1 
0 = 1 
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dp dp dp dp { du dv dw\ 

dt dx dy dz \ dx dy dz j 


dV 


du ( du du 

dw2 ( du 
dy \dy 
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^ dw2 dw2 du 


dx 


dx dx 
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where 


^xz^z^ ^y — ^xy^x ^y^y ^ ^zy^z 


dT dT dT 

^z — ^xz^x ~l“ ^yz^y ^“ ^z^zj 9 ” ^z ^x “I" ^y H" ^z~^ ^z 

ax ay az 


FINITE ELEMENT FORMULATION 


n n 

Let p = Y^ pjxl;j{x,y,z), u = ^ Ujijj{x,y,z), etc. 
j=i y=i 


Equations (1) - (5) can be formulated as 


[M^m + [K^\{p} = {F^} 
[M^]{U} + [K^]{U} = {F^} 
[M^]{V} + [K^]{V}^{F^} 
[M^]{W} + [K'^]{W} = {F*} 
[M^]{t} + \K^]{T} = {F^} 


where 


Mi = K}, = - u,||) ciV 


-L 


, . du dv dw . 

q€ \ ax ay oz 
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= / 


pxliiisjdV, Kfj = 


I / / 



dll)] dtiij 
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dipt d'tbj 
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dy dy dz dz dx dx 
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f \ ^ 
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ALTERNATIVE (CONSERVATION) FORM OF EQUATIONS 


Let V = pv {U = pu, V = pv, W = pw) 
E = pe, /= 0, Q = 0 


The governing equations are 


dp d . d, d . . 


dU 


dp d(Tj 




da 
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dt dx 
dV 
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dt dx'' ' dy 
The finite-element equations are 


lM|W + [iflW = {f‘} 

iMi{(/} + |ir|{(7> = {fn 

[M|{VV} + (if|{»'} = {fn 

iiVf!{£:} + [/f|{£} = {f®} 
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where 


M.y = 


d d d 


[ ti^ii^jdV, Kij = f V't 
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? = / 
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5y 


dV 


dV 


This formulation is a natural extension of the finite element model for inviscid flows and is 
applicable for compressible viscous flows from low subsonic to supersonic flows with suitable 
addition of stabilizing terms (artificial viscosity). For highly viscous, low Mach number 
internal flows there is no need for the addition of artificial viscosity. . This formulation is 
coded in the computer program COMPR3D and is listed in Appendix II. 
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